{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nFollow particle\n===============\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import re\n\nfrom matplotlib import colors\nfrom matplotlib import pyplot as plt\nfrom matplotlib.animation import FuncAnimation\nfrom numpy import arange, meshgrid, ones, unique, zeros\n\nfrom py_eddy_tracker import start_logger\nfrom py_eddy_tracker.appli.gui import Anim\nfrom py_eddy_tracker.data import get_demo_path\nfrom py_eddy_tracker.dataset.grid import GridCollection\nfrom py_eddy_tracker.observations.groups import particle_candidate\nfrom py_eddy_tracker.observations.network import NetworkObservations\n\nstart_logger().setLevel(\"ERROR\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "class VideoAnimation(FuncAnimation):\n    def _repr_html_(self, *args, **kwargs):\n        \"\"\"To get video in html and have a player\"\"\"\n        content = self.to_html5_video()\n        return re.sub(\n            r'width=\"[0-9]*\"\\sheight=\"[0-9]*\"', 'width=\"100%\" height=\"100%\"', content\n        )\n\n    def save(self, *args, **kwargs):\n        if args[0].endswith(\"gif\"):\n            # In this case gif is used to create thumbnail which is not used but consume same time than video\n            # So we create an empty file, to save time\n            with open(args[0], \"w\") as _:\n                pass\n            return\n        return super().save(*args, **kwargs)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n = NetworkObservations.load_file(get_demo_path(\"network_med.nc\")).network(651)\nn = n.extract_with_mask((n.time >= 20180) * (n.time <= 20269))\nn = n.remove_dead_end(nobs=0, ndays=10)\nn.numbering_segment()\nc = GridCollection.from_netcdf_cube(\n    get_demo_path(\"dt_med_allsat_phy_l4_2005T2.nc\"),\n    \"longitude\",\n    \"latitude\",\n    \"time\",\n    heigth=\"adt\",\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Schema\n------\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(12, 6))\nax = fig.add_axes([0.05, 0.05, 0.9, 0.9])\n_ = n.display_timeline(ax, field=\"longitude\", marker=\"+\", lw=2, markersize=5)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Animation\n---------\nParticle settings\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t_snapshot = 20200\nstep = 1 / 50.0\nx, y = meshgrid(arange(20, 36, step), arange(30, 46, step))\nN = 6\nx_f, y_f = x[::N, ::N].copy(), y[::N, ::N].copy()\nx, y = x.reshape(-1), y.reshape(-1)\nx_f, y_f = x_f.reshape(-1), y_f.reshape(-1)\nn_ = n.extract_with_mask(n.time == t_snapshot)\nindex = n_.contains(x, y, intern=True)\nm = index != -1\nindex = n_.segment[index[m]]\nindex_ = unique(index)\nx, y = x[m], y[m]\nm = ~n_.inside(x_f, y_f, intern=True)\nx_f, y_f = x_f[m], y_f[m]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Animation\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cmap = colors.ListedColormap(list(n.COLORS), name=\"from_list\", N=n.segment.max() + 1)\na = Anim(\n    n,\n    intern=False,\n    figsize=(12, 6),\n    nb_step=1,\n    dpi=60,\n    field_color=\"segment\",\n    field_txt=\"segment\",\n    cmap=cmap,\n)\na.fig.suptitle(\"\"), a.ax.set_xlim(24, 36), a.ax.set_ylim(30, 36)\na.txt.set_position((25, 31))\n\nstep = 0.25\nkw_p = dict(nb_step=2, time_step=86400 * step * 0.5, t_init=t_snapshot - 2 * step)\n\nmappables = dict()\nparticules = c.advect(x, y, \"u\", \"v\", **kw_p)\nfilament = c.filament(x_f, y_f, \"u\", \"v\", **kw_p, filament_size=3)\nkw = dict(ls=\"\", marker=\".\", markersize=0.25)\nfor k in index_:\n    m = k == index\n    mappables[k] = a.ax.plot([], [], color=cmap(k), **kw)[0]\nm_filament = a.ax.plot([], [], lw=0.25, color=\"gray\")[0]\n\n\ndef update(frame):\n    tt, xt, yt = particules.__next__()\n    for k, mappable in mappables.items():\n        m = index == k\n        mappable.set_data(xt[m], yt[m])\n    tt, xt, yt = filament.__next__()\n    m_filament.set_data(xt, yt)\n    if frame % 1 == 0:\n        a.func_animation(frame)\n\n\nani = VideoAnimation(a.fig, update, frames=arange(20200, 20269, step), interval=200)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Particle advection\n^^^^^^^^^^^^^^^^^^\nAdvection from speed contour to speed contour (default)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "step = 1 / 60.0\n\nt_start, t_end = int(n.period[0]), int(n.period[1])\ndt = 14\n\nshape = (n.obs.size, 2)\n# Forward run\ni_target_f, pct_target_f = -ones(shape, dtype=\"i4\"), zeros(shape, dtype=\"i1\")\nfor t in arange(t_start, t_end - dt):\n    particle_candidate(c, n, step, t, i_target_f, pct_target_f, n_days=dt)\n\n# Backward run\ni_target_b, pct_target_b = -ones(shape, dtype=\"i4\"), zeros(shape, dtype=\"i1\")\nfor t in arange(t_start + dt, t_end):\n    particle_candidate(c, n, step, t, i_target_b, pct_target_b, n_days=-dt)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(10, 10))\nax_1st_b = fig.add_axes([0.05, 0.52, 0.45, 0.45])\nax_2nd_b = fig.add_axes([0.05, 0.05, 0.45, 0.45])\nax_1st_f = fig.add_axes([0.52, 0.52, 0.45, 0.45])\nax_2nd_f = fig.add_axes([0.52, 0.05, 0.45, 0.45])\nax_1st_b.set_title(\"Backward advection for each time step\")\nax_1st_f.set_title(\"Forward advection for each time step\")\nax_1st_b.set_ylabel(\"Color -> First target\\nLatitude\")\nax_2nd_b.set_ylabel(\"Color -> Secondary target\\nLatitude\")\nax_2nd_b.set_xlabel(\"Julian days\"), ax_2nd_f.set_xlabel(\"Julian days\")\nax_1st_f.set_yticks([]), ax_2nd_f.set_yticks([])\nax_1st_f.set_xticks([]), ax_1st_b.set_xticks([])\n\n\ndef color_alpha(target, pct, vmin=5, vmax=80):\n    color = cmap(n.segment[target])\n    # We will hide under 5 % and from 80% to 100 % it will be 1\n    alpha = (pct - vmin) / (vmax - vmin)\n    alpha[alpha < 0] = 0\n    alpha[alpha > 1] = 1\n    color[:, 3] = alpha\n    return color\n\n\nkw = dict(\n    name=None, yfield=\"longitude\", event=False, zorder=-100, s=(n.speed_area / 20e6)\n)\nn.scatter_timeline(ax_1st_b, c=color_alpha(i_target_b.T[0], pct_target_b.T[0]), **kw)\nn.scatter_timeline(ax_2nd_b, c=color_alpha(i_target_b.T[1], pct_target_b.T[1]), **kw)\nn.scatter_timeline(ax_1st_f, c=color_alpha(i_target_f.T[0], pct_target_f.T[0]), **kw)\nn.scatter_timeline(ax_2nd_f, c=color_alpha(i_target_f.T[1], pct_target_f.T[1]), **kw)\nfor ax in (ax_1st_b, ax_2nd_b, ax_1st_f, ax_2nd_f):\n    n.display_timeline(ax, field=\"longitude\", marker=\"+\", lw=2, markersize=5)\n    ax.grid()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}